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GLOBAL EQUATION SOLVER AND OPTIMIZER 
Background 

[0001] Planning and scheduling require solving and optimizing a set of nonlinear 

equations. Consider a plant that makes ice cream. The plant needs to make 
vanilla, strawberry, and chocolate ice cream, but there are a number of 
constraints, such as the amount of vanilla flavoring, strawberries, and chocolate 
available. There are scheduled deliveries of ice cream and orders to fill. When a 
truck rolls up and expects to be loaded with chocolate ice cream, for example. 
There are production constraints, such as needing to do chocolate last because 
after producing chocolate ice cream, the pipes must be cleaned before you can 
run anything else. This is a cost penalty for chocolate. While the plant is 
making one kind of ice cream it cannot produce anything else, because the 
equipment is committed to a particular product. Business managers would like 
to make as large a batch as possible because that minimizes the amount of time 
the system is down for transitions between products. On the other hand, storage 
space is limited. For example, chocolate syrup may start to build up. Trucks full 
of chocolate syrup sit in the plant parking lot and charge thousands of dollars a 
day until you can get them unloaded. 
[0002] In managing plant operations, there are a series of continuous and 

discreet decisions to make. To get rid of the chocolate syrup, chocolate will need 
to be produced on Tuesday and strawberry on Wednesday and some time will 
need to be spent cleaning up after the chocolate. Those are discreet decisions, 
because a plant cannot only half make chocolate ice cream; the plant either 
makes it or not. Given a particular set of discrete decisions, a set of continuous 
decisions is made. How much chocolate ice cream should be produced? How 
much strawberry ice cream should be produced? How much chocolate is the 
plant going to have in the tanks at a particular time of day? How much chocolate 
will be produced in the next hour? Is it going to build up to the point where is 
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blows the top off the tank? How much chocolate syrup is required? When does 
the plant start making strawberry? How fast will the strawberries be used up and 
when will the strawberries run out? Those are continuous decisions. 
Complications may arise, such when the guy who makes strawberry is late or 
when it takes longer than usual to clean up after the chocolate. Discrete and 
continuous decisions like these are used to set up a set of equations or 
constraints. 

[0003] In addition, there are other motivations, such as market demand. A 

customer may be willing to pay twice as much for all chocolate ice cream as for 
half chocolate and half strawberry. Perhaps a plant manager is motivated to do a 
lot more vanilla, even though the supply schedule tells him he is not able to get 
to vanilla for a long time. Then, he will want to optimize production of vanilla 
as soon as the vanilla supplies arrive. Even more complex situations occur. 

[0004] For a planning system, a set of equations derived from high-level 

constraints is solved to find the best solution based on certain criteria, such as 
profitability. This result tells whether the plan will work. If so, then more 
detailed constraints and criteria are added to the set of equations and the set of 
equations is solved for a schedule. Scheduling systems have a much larger set 
of equations to solve than planning systems. This puts a lot more stress on the 
solver and the solving technology. 
[0005] Most current systems are planning systems, not scheduling systems, 

which are larger and harder to solve. Current methods are unable to solve 
scheduling problems within the necessary time for rapid modeling and 
simulation with optimum or near optimum solutions. Also, current methods 
suffer from a local optimality problem, where a local optimum in the solution 
space is incorrectly given as the solution when the actual global optimum lies 
elsewhere. Any system that suffers from the local optimality problem may not 
be able to find a solution, when a solution exists. In addition, some of the 
current methods are too slow. Furthermore, most current methods fail to find a 
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solution without indicating if the problem is infeasible or not. The few methods 
that are able determine infeasibility have poor or slow convergence. Some 
examples of current methods are sequential linear programming, sequential 
quadratic programming, and trusted regions. All three of these suffer from the 
local optimality problem and when they fail, it is unknown whether the problem 
is infeasible or not. There is a need for an efficient method to quickly converge 
on a solution, to rigorously determine whether or not the problem is infeasible, 
and to find the global optimum instead of getting stuck in a local optimum. 

Summary 

[0006] The present invention solves both planning and scheduling problems by 

combining a number of technologies. It balances the safety of subdivision 
methods with the fast convergence of linearization methods, avoiding the local 
optimality problem. Also, the linearization methods rigorously determine 
whether or not the problem was infeasible. The present invention is a method of 
solving an operations problem. Operations problems comprehend both planning 
and scheduling problems. First, the method receives variables, relationships, and 
constraints. Then, it forms a set of non-convex quadratic equations based on 
them. The method solves the set of non-convex quadratic equations by applying 
a bound propagation process, a local linear bounding process, a local 
linearization process, and a global subdivision search. Finally, the method 
determines whether a solution is optimal, feasible, or infeasible. Thus, the 
present invention recognizes local optimality problems and goes beyond them to 
find a global solution, if one exists. If not, the present invention rigorously 
proves infeasibility. If there is a solution, it comes up with a solution. 

Brief Description of the Drawing s 
[0007] Figure 1 is a block diagram of example applications of the present 

invention. 
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Figure 2 is a block diagram of a bounded region containing a solution 
space to solve for an optima using embodiments of the present invention. 

Figure 3 is a flow chart of a method embodiment of the present invention. 

Figure 4 is a more detailed flow chart than Figure 3 and shows a method 
embodiment of the present invention. 

Figure 5 is a block diagram of a bounded region, like the one shown in 
Figure 2, which is modified in a bound propagation and refinement subprocess of 
a method embodiment of the present invention, such as the method embodiments 
shown in Figures 3 and 4. 

Figure 6 is a flow chart of a local linear bounding subprocess of a method 
embodiment of the present invention, such as the method embodiments shown in 
Figures 3 and 4. 

Figure 7 is a flow chart of a linearization subprocess of a method 
embodiment of the present invention, such as the method embodiments shown in 
Figures 3 and 4. 

Detailed Description 

[0008] Methods for a global equation solver and optimizer are described. In the 

following detailed description, reference is made to the accompanying drawings, 
which form a part hereof. These drawings show, by way of illustration, specific 
embodiments in which the invention may be practiced. In the drawings, like 
numerals describe substantially similar components throughout the several 
views. These embodiments are described in sufficient detail to enable those 
skilled in the art to practice the invention. Other embodiments may be utilized 
and structural, logical, and electrical changes maybe made without departing 
from the scope of the present invention. 

[0009] Figure 1 is a block diagram of example applications of the present 

invention. The present invention is shown as a solver 100. The solver solves a 
set of equations, giving a set of variables and what they are equal to in the 
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solution, if one exists. These variables are then translated into a business use, 
such as an operations plan, an accounting summary or some other useful, 
concrete, and tangible results. The solver is like an octopus, because customer 
orders interact with inventory, which in turn interacts with production facilities, 
which interacts with the business plan and so on. The solver is at the heart of all 
these activities. 

[0010] The solver is run as a scheduling system at least on a daily basis. As a 

planning system, it is run about once a month with adjustments as events occur. 
Or run to interact with traders as they phone in opportunities. The solver 100 
can be applied to scheduling 102, planning 104, and other operations 106. 

[001 1] Planning is the forecasting of the average performance of a plant (a 

collection of interconnected processes) over some specified period, such as a 
month or a year. The plan specifies what inputs are needed and how they are to 
be used to produce plant outputs. The plan usually includes forecasts of the 
values of individual process performance parameters such as yields, product 
qualities, flow rates, temperatures, and pressures. The plan also includes 
economic information about the impact of changes in parameters on plant 
profitability. 

[0012] Scheduling is the specification of the inputs to and outputs from each 

process and inventory, plus the timing and sequencing of each production 
operation, whether batch or continuous, over some short scheduling period, such 
as a week or 10 days. Although the horizon is a week or 10 days, today's 
operation is the most important. Operations are not averaged over the scheduling 
period; rather, time and operations move continuously from the beginning of the 
period to the end. Ideally, the schedule is revised each day as needed so that it 
always starts from current actual performance. 

[0013] Additionally, the solver can be applied to inventory 108, suppliers 1 10, 

and ordering 1 12. Also, it can be applied to customers 1 14 and production 116. 
For example, a business analyst takes a phone call from a trader who says he has 
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a chance to sell this much of something and asks the business analyst if he can 
make it in time. The business analyst quickly runs the solver to decide whether it 
is a good opportunity or not. Another example is storing routine results in an 
executive information system for general business planning, such as forecasting 
whether profit numbers will be made for the quarter. Another example is an 
operations person who uses the results to decide what temperature to crank his 
unit. Another example is a factory floor scheduler with a cell phone and 
spreadsheet, who directs operations people based on the results. In addition, the 
present invention has many other applications. 
[00 1 4] Figure 2 is a block diagram of a bounded region 200 containing a 

solution space 202 to solve for an optima 204 using embodiments of the present 
invention. Initially, the solution space is unknown, but the bounded region 200 
can be determined from the set of equations to solve. The ranges of variables in 
the set of equations define certain bounds. The exemplary bounded region 200 is 
shown in three-dimensional space with x- 206, y- 208, and z-axes 210. This 
represents a rather simple set of equations with 3 variables. A typical set of 
equations has 10,000 variables, so a typical bounded region would have 10,000 
dimensions. However, the three dimensions shown in Figure 2 are easier to 
conceptualize. In Figure 2, the x variable has lower and upper bounds, LB(X) 
and UB(X) respectively. The y variable has lower and upper bounds, LB(Y) and 
UB(Y) respectively. The z variable has lower and upper bounds, LB(Z) and 
UB(Z) respectively. The bounded region 200 is the intersection of these lower 
and upper bounds and the solution space 202 is contained within the bounded 
region 200. 

[00 1 5] Assuming the problem is not infeasible, the solution space 202 is a set of 

feasible solutions in a feasible region, one of which is the optima 204. A 
feasible region is the set of all points satisfying all of the problem's constraints 
and restrictions defined in the set of equations. For example, in a feasible 
solution for an oil refinery, no tanks overflow or underflow. Often, for a 
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scheduling problem, the goal is to meet the schedule with a feasible solution, 
rather than an optimal solution. For example, given a certain amount of 
inventory and shipments, a feasible solution makes the shipments without using 
up too much of the inventory. Because the deals have already been made, any 
extra product is stored, not sold. An infeasible problem occurs when the feasible 
region is empty, i.e. it contains no points. By definition, an infeasible problem 
has no optimal solution. An optima or optimal solution is a point in the feasible 
region with the largest objective function value, for a maximization problem. 
Similarly, for a minimization problem, the optimal solution is a point in the 
feasible region with the smallest objective function value. An objective function 
is some function of particular decision variables to be minimized or maximized. 
For example, a simple objective function to maximize is (weekly revenues) - 
(raw material purchase costs) - (other variable costs). Some other examples of 
optimality measures or objective functions are minimum turnover, minimum 
costs, accomplishing work in minimum time, maximizing a high margin product, 
and the like. The optima is defined by pre-determined criteria, such as least cost 
or most profit that are included as input to the present invention. A global 
optima is the optima over the entire range of variables as opposed to a local 
optima, which is the optima only in a local area. 
[001 6] Figure 3 is a flow chart of a method embodiment 300 of the present 

invention. One aspect of the present invention is a method 300 of solving an 
operations problem. Operations problems comprehend both planning and 
scheduling problems. For example, problems include maximizing profits, 
meeting shipments, meeting production schedules, meeting product specification 
requirements, and other business problems. Operations problems are 
optimization problems in industries as diverse as banking, education, forestry, 
petroleum, and trucking. The method 300 comprises receiving variables 302, 
relationships 304, and constraints 306. The variables 302 are things like 
qualities, quantities, timing, and the like. Relationships 304 express how the 
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variables 302 interrelate, such as how speed relates to quality. For example, in 
mixing things, the relative quantities effect the quality of the mixture. Other 
examples are relationships 304 based on the physics of a plant and relationships 
304 based on economics, such as cost and revenue. Some examples of 
constraints 306 in refinery applications are tank limits, product specifications, 
gasoline octane ratings, operating limits and the like. 

[001 7] The method 300 comprises forming a set of non-convex quadratic 

equations 308 based on the variables 302, relationships 304, and constraints 306. 
Some example equations are formed by the definition of the problem. Other 
example equations are formed by physical constraints, such as the capacity of a 
tank. Still other example equations are formed by the physics of a process being 
modeled. Equations that are not in quadratic form can be converted to quadratic 
form by standard approximation methods. Convex equations have solution sets 
such that a line segment joining any pair of points in the solution set is wholly 
contained in the solution set. Non-convex equations do not have this property. 
Generally, when a convex object, such as a beach ball is put on a flat surface, it 
will be touching at one point. Non-convex objects basically have some 
indentation or dimple in it. Non-convex equations usually have multiple local 
pseudo-solutions, meaning they are not solvable with current methods having the 
local optimality problem. 

[001 8] The method 300 further comprises solving the set of non-convex 

quadratic equations by applying a bound propagation process, a local linear 
bounding process, a local linearization process, and a global subdivision search. 
This is shown in Figure 3 as the solver 310. Solving the equations is described 
in more detail below with reference to Figure 4. The method 300 further 
comprises determining whether a solution is optimal, feasible, or infeasible. 
Figure 3 shows the solver 310 determines whether a solution exits 3 12. If no 
solution exists, then the solver 310 determines that the problem is infeasible 314. 
If a solution exists, then the solver 310 determines whether the solution is 
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optimal 3 1 6. The solver determines either that the solution is optimal 3 1 8 or 
feasible 320. 

[001 9] The solution to the set of non-convex quadratic equations has many uses. 

(See Figure 1 .) In one embodiment, the solution is a schedule for the 
manufacturing process. In another embodiment, the solution is a schedule for 
operating an oil refinery. In another embodiment, the solution is a plan for the 
manufacturing process. In another embodiment, the solution is a plan for 
operating an oil refinery. 
[0020] Figure 4 is a more detailed flow chart than Figure 3 and shows a method 

embodiment 400 of the present invention. One aspect of the present invention is 
a method 400 comprising solving the set of non-convex quadratic equations by 
applying a bound propagation process 402, a local linear bounding process 404, a 
local linearization process 406, and a global subdivision search 408. Three of 
these processes are described below with reference to Figures 5-7. The bound 
propagation process 402 is described below with reference to Figure 5. The local 
linear bounding process 404 is described below with reference to Figure 6. The 
local linearization process 406 is described below with reference to Figure 7. 
[002 1 ] The global subdivision search 408 is done over a bounded region, like the 

one shown in Figure 2 using branch and prune logic. For branching, the solver 
starts with the whole bounded region and if no solution is found, splits the region 
in half, first searching one half and then the other until a solution is found or the 
problem is determined to be infeasible. For pruning, only bounded regions that 
might contain a solution need to be examined. There are many alternate methods 
for performing the global subdivision search. One alternate method is to pick the 
most infeasible constraint and find the most infeasible variable in it. 
[0022] The combination of processes in method 400 work together to create 

advantages over current methods. The bound propagation process 402 gives the 
method 400 improved performance over current methods. The local linear 
bounding process 404 rigorously proves infeasibility locally. The local 
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linearization process 406 improves time to convergence on a solution over 
current methods. The global subdivision search 408 helps to avoid the local 
optimality problem. These advantages help reduce project cycles, integrate 
operations for better decisions, and increase profits by providing more accurate, 
timely information for real-time decisions. 

[0023] Another aspect of the present invention is a machine-accessible medium 

having associated content capable of directing the machine to perform a method 
400 of solving a set of non-convex quadratic equations. At the start 41 0 in 
Figure 4, the method 400 comprises selecting a region bounding all variables 
412. A bound propagation process is applied to the region to refine the bounds 
and improve linearization 402. A local linear bounding process is applied to the 
region to determine feasibility and to find approximately feasible solutions 404. 
The local linear bounding process is optionally repeated. The method 400 
determines if there is no potential for a solution 414 within the region. If so, the 
region is eliminated from consideration 416. The method 400 determines if 
there is no potential for an optimal solution 41 8. If so, the region is eliminated 
from consideration 420. 

[0024] A local linearization process is applied to the region to determine 

feasibility and local optimality 406. The local linearization process is optionally 
repeated. The method 400 determines if there is a solution 422 and if so, if it is 
optimal 424. Upon finding an optimal global solution 426, the optimal global 
solution and information indicating optimality are provided and the method 400 
ends 428. Upon finding a feasible global solution 430, the feasible global 
solution and information indicating feasibility are provided and the method 400 
ends 428. Upon determining local infeasibility, the region is eliminated from 
consideration. Upon determining global infeasibility 432, information indicating 
infeasibility is provided and the method 400 ends 428. Upon not finding a 
solution, a global subdivision search 408 is applied to the region to produce two 
or more regions and the bound propagation 402, local linear bounding 404, and 
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local linearization 406 processes are iteratively applied to each of the two or 
more regions, until the solution is determined to be optimal 426, feasible 430, or 
infeasible 432. 

[0025] In one embodiment, the method 400 further comprises receiving input 

variables, constraints, and equations. In another embodiment, the method 400 
further comprises receiving a measure of optimality used to determine the global 
optimal solution. In another embodiment, the method 400 further comprises 
receiving a measure of feasibility used to determine the global feasible solution. 
In another embodiment, the method 400 further comprises providing a schedule 
for operating a plant. In another embodiment, the method 400 further comprises 
providing a plan for operating a plant. 

[0026] Another aspect of the present invention is a process 400 of solving a set 

of non-convex quadratic equations. The process 400 comprises selecting a 
region bounding all variables 412. A bound propagation process 402 is applied 
to the region to refine the bounds and improve linearization. A local linear 
bounding process 404 is applied to the region to determine feasibility and to find 
approximately feasible solutions. A local linearization process 406 is applied to 
the region to determine feasibility and local optimality. Upon finding a solution 
after the local linearization process, the solution is provided. Upon determining 
infeasibility, the region is eliminated from consideration. Upon not finding the 
solution after the local linearization process, a global subdivision search 408 is 
applied to the region to produce two or more regions and the bound propagation 
402, local linear bounding 404, and local linearization 406 processes are 
iteratively applied to each of the two or more regions, until it is determined that 
the solution is optimal, feasible, or infeasible. 

[0027] In one embodiment of the process 400, the local linearization process 406 

is the local linear bounding process 404. 

[0028] Figure 5 is a block diagram of a bounded region, like the one shown in 

Figure 2, which is modified in a bound propagation and refinement subprocess of 
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a method embodiment of the present invention, such as the method embodiments 
shown in Figures 3 and 4. Figure 5 shows a bounded region 500 with refined 
upper and lower bounds along the y-axis to produce a refined region 502, which 
is closer to the solution space 504 and reduces the space to search. Once bounds 
are refined, they are propagated. For example, if one of the equations is X = Z + 
Y, once Y is bound, then that can be propagated to X and Z so that X and Z are 
also bound, reducing the size of the region to search. 
[0029] Figure 6 is a flow chart of a local linear bounding subprocess 600 of a 

method embodiment of the present invention, such as the method embodiments 
300, 400 shown in Figures 3 and 4. In one embodiment, the local linear 
bounding subprocess 600 determines infeasibility and finds approximately 
feasible solutions. In one embodiment, the local linear bounding process 600 
comprises performing differentiation 602 on equations 604 in the region 606. 
The region 606 includes an initial point (x 0 ). Differentiation 602 is done around 
the x 0 . The lower and upper bounds on the variables in the region are 
determined 608. A linear programming process is applied to the linear equations 
in the region 610. The local linear bounding subprocess 600 determines whether 
a solution exists in the region 612. Upon finding a solution exists, local 
feasibility is determined 614. Upon finding local infeasibility, global 
infeasibility is determined 616. 
[0030] In one embodiment, the initial point is a trial solution closest to a feasible 

solution. In another embodiment, the initial point is a trial solution closest to 
being the optimal solution. In another embodiment, the initial point is the 
largest. In another embodiment, the initial point is the smallest. Other 
embodiments include using discrete search techniques. 
[003 1] Figure 7 is a flow chart of a linearization subprocess 700 of a method 

embodiment of the present invention, such as the method embodiments 300, 400 
shown in Figures 3 and 4. In one embodiment, the local linearization process 
700 comprises performing differentiation 702 at a point in the bounded region 
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704. A set of linear equations is formed 706. A linear programming process is 
applied to the linear equations in the bounded region 708. A new point is 
generated in the bounded region 710 and the local linearization process is 
repeated with the new point. If the linear program fails to generate a new point 
712, then the global subdivision search or some other method is applied 714. In 
another embodiment, the local linearization process 700 has quadratic 
convergence, giving it fast performance. 



solution and optimization of systems of quadratic equations, with extensions to 
general nonlinear functions. 



[0032] 



Example Embodiments 
The following example embodiments provide methods for the global 



Basic Definitions 



[0033] 



Let f(x) : 9T 9T be a quadratic function of the form 



(1.1) AW^Ct+Z^Xi+X^XjXi- 



The problem we wish to solve will have the following fonn. 



min f k (x):keK 0 



(1 .2) lb k < f k (x) < ub k : Vk e K x 



(The set K 0 must include only a single linear function.) 



[0034] 



With informal notation, we shall write the functions in the matrix/vector 
notation 



(1.3) f(x) = C + Ax + Bxx 
when we wish to consider all the functions, and in the forms 




[0035] 



f l (x) = C l + A x x + B l xx 
when we wish to consider the three classes of functions. 

Notationally, Bxy is to be interpreted as a multilinear 
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form Bxy = J] B fg . i x j y i and the linear form Bx is equivalent to the matrix 

y 

== Yj B ^ x j ' 131 Edition, the transpose of B is defined as B *xy = Byx , 

y 

with the immediate equality B *^ = 5 fcy . 

[0036] Define the positive and negative functions 

(15) = max(x,0) >0 

(*)_ = min(x,0) < 0 

[003 7] We define a measure of the infeasibility of a particular point to be 

(1.6) AW=jA t (x) 

where 

(1.7) A, (*) = max((/ ft (*) - ub k ) + , (ft 4 - /, (,)) + ) . 

(As a result, A(x) > 0 , and A(x) = 0 only if the point ;c is feasible.) 
[0038] In the course of solving the problem (1 .2) above, we will be solving a 

sequence of subsidiary problems. These problems will be parameterized by a 
trial solution and a set of point bounds, 

(1.7) * 

u<x<v 

From that we can compute a set of gradient bounds 

(18) F = W+ U + ( B )- V 
G = (B) + v + (B)_u 

so that 

(1.9) u<x<v => F<Bx<G. 

We will also require that the trial solution also satisfy the bounds 

(1.10) w<3c<Cv. 

[0039] We define the row vector e = {1,...,1} so that we can compute the 

maximum infeasibility for the bounds 

(1.11) A(a,v) = e(G-F){v-u) = £(G tt -F u )(v. -w,.) 
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(1.12) A(«,v) = e \B\(v-u){v-u) = ^JfVj - Uj ){ Vi -u ( ) . 

kii 

[0040] In the course of the method, we will also be computing a series of trial 

solutions. We can compute the divergence of a sequence of trial solutions 

d.13) i^-^i^a^-^i. 

i 

The centered representation of a function relative to a given trial solution 3c is 
(1.14) f( x ) = C + A(x-x) + B(x-x)(x-x) 

where 



A = A + (B + B*)x 

(1.15) 



V 



f 

C =C + Ax + Bxx 



V 



[0041] By also defining 



u =u-x 

v - v-x 

(1.16) _ 
F=F-Bx 

G=G-Bx 

the bounding inequalities are equivalent to the following centered inequalities 

u < x-x <v 

(1.17) _ 
F<B(x~x)<G 

[0042] Also note that since u < 0 < v , we have 

(1.18) |*. - x.\ < maxd^.yvl.) < v, -«. . 

In order to bound on both side, we will be decomposing x-x into two 
nonnegative variables, 

z-w=x~x 
(1.18) z>0 
w>0 
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to which we will apply the bounds 

(1.19) z.+w, <max(|w,|,|v|.). 

[0043] The subsidiary problems are then 

A) the basic quadratic feasibility problem- 



(1.20) P0 = \x: 



u°<x<v° 



lb<f\x)<ub\ 



and its optimization form 

(1.21) POO J, : min /°Wl 
[ xePO j 

B) the basic quadratic problem with the bounding inequalities- 

u < x < v ] 



(1.23) Pbd{u,v) = \x: 
and its optimization form 

(1.24) PObd(u,v)^\x: 



lb<f\x)<ub\ 



mmf\x) 



xe Pbd{u,v)\ 

C) the enveloping linear programming problem (using the formulas of (1 .8), 
(1.15), and (1.16)) — 



(1-25) LP(x,u,v) = < 



x,z,w: 



U < X ~ X < V 

lb<C l +A\x-x) + G l z-F l w 

ub>C l + A\x-x) + F l z-G l w 

x-x - z-w 

z + w<max(|vj,|ttj) 

z,w>0 



its optimization form 

(1.26) LPO(x,u,v) = ^x: 
and its minimal infeasibility form 

(1.27) LPMI(x,u,v) = \x: 



mmf°(x) = A°x} 
xeLP(x,u,v) I 



mine(G-F)(z+ w)l 
xeLP(x,u,v) I 
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E) and finally the linearization problem with the bounding inequalities included 

min A 0 x 

x: u <x-x <v 

lb<C l +A\x-x)<ub 



(1.28) LLP(x,u,v) = 



Method 

[0044] The problems above have the following relationships. 

Pbd(u 9 v)QLP(x,u,v) 
If x r eLP(x,u,v) then 

Pbd(u, v) n LP{ x, u, v) = Pbd(u 9 v) n LP(x\ u 9 v) 
If x f e ZP(x,a,v)then 

A(x') < e(G - F)(z' + W) < A(u, v) = e(G - F)(v - u) 

(Note that it is also true that if G H -F H = 0 , then G w = F h = {Bx) H = 0 , and 

that particular quadratic term will disappear from the constraint approximations 
in LP(x,u,v). ) 

[0045] As we search for a solution for a problem POO , we split the region up 

into nodes, each of which is given by a set of bounds {u 9 v} . 

[0046] For any problem, the theory of Newton's method tells us that there is a 

constant 0 > 0 such that for any A(w, v) < 9 , if Pbd(u 9 v) is feasible then the 
sequence of trial solutions generated by the linearization, x m+l = LLP(x m 9 u,v) , 
will converge quadratically to the solution ** = Pbd(u 9 v) , with 

[0047] On the other hand by the theorems proven here, there is a constant 

<9 > 0 such that for any A(w, v) < 0 , if Pbd(u,v) is infeasible then LP(x,u 9 v) 
will be infeasible. 

[0048] So a simplistic view of the method is to 

1) Choose a node {u 9 v) , and find a trial solution u < x < v , 

2) Run J"" 1 = LLP(x m 9 u 9 v) and see if it finds a feasible point, 
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3) Run x € LP(x,u,v) and see if it is infeasible, 

4) If the node fails both (2) and (3), subdivide the node into smaller regions, and 
try again. 

Expanding a Node 

[0049] An open node specifies a particular trial solution and bounds, {x,u y v} . It 

also specifies a particular distribution of the quadratic terms B + B * (see the 
Symmetry Breaking of the Gradient Terms section below). 

[0050] To expand the node, apply the procedure below until the node is fully 

expanded, or has been closed. The procedure will update the problem, its trial 
solution and bounds, will establish the status of the node, and will update a 
rigorous lower bound cp of the objective function on the node. 

[005 1 ] Note that in the course of searching for a suitable node to expand (see 

section below), it may be desirable to apply only steps (1), (2), (3), and (4) in 
order to see if there exists a linearization node which succeeds in step (2). This 
will be referred to as partial expansion. 

[0052] There are two basic patterns of application of these steps — 

A) Run linearization alone as long as possible — 

Run steps (1), (2), (3), and (4) when partially expanding a node. Run 
steps (5) and (6) when fully expanding a node. 

B) Always seed linearization with a minimal infeasibility trial solution- 

Run steps (1), (4), (5), (2) and (3) when partially expanding a node. Run 
steps (5) (rerun) and (6) when fully expanding a node. Step (2) 
should include multiple iterations of linearization, otherwise the 
quadratic convergence property will be lost. 
[0053] Let s > 0 be the basic numerical tolerance desired. 

[0054] 1) Propagate the bounds through the problem (see the Convergence and 

Divergence of Trial Solutions section below). Update the node with the new 
bounds and/or problem terms — 
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{u,v} = 

{C,A,B} = {CU',B'y 
[0055] 2) Try to solve the linearization problem LLP(x, w, v) . 

Optionally, one can solve a series of linearization problems, 
x m+l = LLP(x m ,u,v) . The series should be ended when: 

The problem becomes infeasible; 

A fixed upper limit of iterations is reached; 

The trial solutions converge; or 

The trial solutions diverge. 
At the end of this series, if there exists a solution jc' e LLP{x,u,v) , then use the 
solution as the new trial solution x = x r , and declare the node linearized. 
[0056] If there were a series of trial solutions found, then if the solutions 

converges or the upper limit reached, choose the final trial solution, but if the 
problem became infeasible or the trial solutions diverged, choose the trial 
solution with the minimum infeasibility, x r = argmin m A(x m ) . 

[0057] Ifthe series was ended because the solutions converged, declare the node 

convergent, 

[0058] Ifthe series was ended because the solutions diverged, declare the node 

divergent. 

[0059] In any case, set the optimal lower bound to the worst bound on the node, 

p = (A°) + u + (A°)_v. (This will be reset in step 6, but is set here to support 
partial expansion.) 

[0060] 3) If step (2) succeeds in determining a trial solution, compute the 

infeasibility of the new trial solution, A(3c) , and the maximum infeasibility of 
the bounds A(w,v). 

[0061 ] If A(w, v) < s , all points satisfying x' e Bd{x, w, v) are feasible within the 

desired tolerance. Declare the node completely feasible and point-optimal. 
If A(x) < s , then the trial solution is feasible within the desired tolerance. 
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Declare the node point-feasible. (Note that completely feasible implies point- 
feasible.) (Note: if the node was convergent, it should be point-feasible.) 
If feasibility is all that is desired (as opposed to optimality), make the following 
substitution for steps (4)-(6>— If the node has been declared point-feasible, then 
close the node. If not, the node is completely expanded. Exit this procedure in 
either case. 

[0062] If the node is linearized and point-feasible, then by the standard theory of 

Newton iteration, the trial solution is an approximate local solution to the 
nonlinear bounded problem Bd(u,v,F y G) . Close the node and declare it point- 
optimal, and set the optimal lower bound for the node to the linearized optimum, 
V=A°x. 

[0063] Alternately, in the case where there are multiple local solutions within the 

same bounds, determine a local region within which the point is optimal, and 
close that newfound node, declaring it as point-optimal and point-feasible. Add a 
constraint to the original node such that the optimum must be strictly better than 
the found local optimum, A°x < cp - s and treat it as a new unopened node. 

[0064] 4) If step (2) fails, try to solve the enveloping minimal infeasibility 

problem LPMI(x,u,v). 

[0065] If no such solution exists, close the node and declare it infeasible. Exit 

this procedure. 

[0066] If there exists a solution x f e LPMI(x,u y v) , then use the solution as the 

new trial solution x-x T . (Also record the auxiliary variables {z\W} .) 
[0067] 5) Compute the infeasibility of the new trial solution, A(x) , and the 

maximum infeasibility of the bounds A(u 9 v) . 
[0068] If A(w, v) < s , all points satisfying x' e LPMI(x,u,v) are feasible within 

the desired tolerance. Declare the node completely feasible. 
[0069] If A(x) < s , then the trial solution is feasible within the desired tolerance. 

Declare the node point-feasible. (Note that completely feasible implies point- 
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feasible.) 

[0070] 6) If the node is not linearized, and is not infeasible, solve the enveloping 

optimal problem, x" e LPO(x,u,v) . Set the optimal lower bound to the resulting 

minimum, (p = A°x n . 
[0071] If the node is completely feasible, then use the new solution as the new 

trial solution x - x" , close the node, and declare it point-optimal. Exit this 

procedure. 

[0072] If the node is point-feasible and the lower bound is close enough to the 

trial solution, J<p - A 0 #| < s , close the node, and declare it point-optimal. (Do 

NOT use the new solution as the trial solution.) Exit this procedure. 
Otherwise, the node is completely expanded. Exit this procedure. 
[0073] Alternately, if feasibility is all that is desired (as opposed to optimality), 

make the following substitution for step (6) — If the node has been declared 
point-feasible, then close the node. If not, the node is completely expanded. Exit 
this procedure in either case. 



Choosing, Expanding, and Splitting a Node 
[0074] Define <p * be the current best optimum found at a particular point in the 

search. That is 

^ = Jmin(^(A0 : Nisa point - optimal node) 
( oo if no point -optimal node exists 

Let s > 0 be the basic numerical tolerance desired. Then we proceed as follows 
with the list of non-closed nodes. 

[0075] 1) For all non-closed nodes N , if cp(N) >q>*-s, close the node and 

declare it sub-optimal. Also, we can close those with bounds equal within 
tolerance to the best as well as those with bounds strictly greater. Alternately, if 
feasibility is all that is desired (as opposed to optimality), this step can be 
skipped. 
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[0076] 2) Choose a candidate set of nodes from the set of non-closed nodes 

according to the following. 

[0077] If there are linearized nodes, select that set and proceed to step (3). 

[0078] If there are nodes that have not yet been partially expanded, then partially 

expand nodes until either one is linearized, or all have been partially expanded, 
then go back to step (1). 

[0079] At this point, one may optionally choose to fully expand all non-closed 

nodes, if there are any that need it, and go back to step (1). 

[0080] If there are no fully expanded non-closed nodes and there are any non- 

closed nodes that have not been fully expanded, expand at least one of them and 
go back to step (1). 

[0081] Select the set of expanded nodes and proceed to step (3). 

[0082] Else all nodes have been closed, and you may terminate the procedure. 

[0083] 3) From the set of non-closed nodes, we wish to select one to split. From 

the set of chosen nodes, select an individual node according to one of the 
following possible measures. 

The node with the smallest infeasibility of the trial solution, A(x) . 
The node with the smallest ^ . 

The node with the largest maximum infeasibility of the bounds, A(w, v) . 

These are suggested heuristics. Correctness only requires the expansion of an 
open node. 

[0084] 4) It is now time to subdivide the node. 

[0085] We will subdivide the point range. Compute w = ( u * v ^ . 
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Optionally, one could also choose w = x if one wants to subdivide at the points 
of linearization. 

[0086] If the node is linearized, and divergent, we will compute the worst 

divergence 

i* = argmax {j . } |x^ +l - 5c* j . 

If the node is not linearized, , we will compute the dimension of the largest 
infeasibility according to the trial solution. 

i* - argmax ( ,. } £(G fa . -F u Xz t + w.). 

k 

Otherwise, we can compute the dimension of the largest infeasibility according 
to the point bounds. 

f* - argmax {/} £(G fe . -F fa .)(v ; -u t ) . 

[0087] 5) Close the node, and open two new nodes with the same problem, trial 

solution, and bounds are the newly closed node, except that for the first new 
open node, 

v' f =v. : Vi */* 

v;» = +s/io 

and for the second new open node, 

u[ = u i : Vi * i * 

w;* = w,*-s/io . 

v;=v,:Vz 

Initialization and Termination 
[0088] To start the list of nodes, we require a problem {C, A, B} , a realistic set 

of finite bounds with a reasonable trial solution - oo<w<3c<v<oo and a basic 
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error tolerance e . 

[0089] To terminate the procedure, we continue to split, expand, and close nodes 

until all nodes have been closed. At that point we have either: 

1) There exists a point-optimal node whose value is as minimal as any other 
node, and thus the trial solution of that node is the optimum solution of the 
original problem; or 

2) All nodes are infeasible and so is the original problem. 

Alternately, if all we are interested in is feasibility, we only need look for any 
node that is point-feasible, in which case the trial solution is the feasible point. 

Convergence and Divergence of Trial Solutions 
[0090] By the theory associated with Newton's method, if the bounded problem 

is feasible, and the initial trial solution is "sufficiently" close to the solution 
x* = Pbd(u,v) 9 then the series of trial solutions will be quadratically convergent, 

with ||x m+2 - < a\\x m+] - x m f . If the problem is infeasible, or if the initial 

trial solution is not sufficiently close, then examples from chaos theory show that 
almost any behavior may occur. 
[0091] So our criterion for convergence or divergence will be the following. 

The trial solutions diverge if \\x m+2 - x m+x || > ||jc W+I - x m | . 

The trial solutions converge if jjjc m+2 - x m+] jj < s . 
Otherwise, we pass no judgment. 

The measures used for convergence/divergence are purely pragmatic, and the 
metric |j| can be chosen for convenience, so that for example either the sum of 
difference, or the maximum difference can be used. 



Simple Propagation of Bounds 
[0092] Iterate the following propagations until no further improvement is seen. 

Alternately, include the square-free bound and the non-square-free bound 
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propagations, until no improvement is found. 
[0093] 1 ) If there exists a {k', i '} : {G M -F^ks, then to within the error 

tolerance, we can assume that 

and we can then remove a quadratic term from the problem by adding the linear 
constraint 



2 



and modifying the original constraint 

f k . (x) = C k . + £ 4,.*. + £%*,* f 

to 

[0094] Let the new constraint have new index k . Then 

C' k =C k :Vk*k' 

* 2 
4 = A u :Vk*k'Vi 

- + 

B' ¥ =B ¥ :Vk*k'VjVi 

%, = 0:Vy 
^ = 0:VyV/ 

[0095] 2) For every linear constraint 
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we can determine new point bounds 



VkeK^ f k linear, Vz : < 



A u > 0 => v[ = min 



H 



4, < 0 => v! = min 



%-c;-E(4) + v;-5:k).«; 

v ;, £ — ^ 



Vfc €K 19 f k linear, V/ : < 



A u > 0 => w' = max 



f %-ci-Z(4) + v;-i(4)-»; 



^ < 0 u z ' = max mJ, 



Symmetry Breaking of the Gradient Terms 
[0096] We can exploit the fact that the original quadratic functions only are 

affected by the value B + B * in order to attempt to reduce the degree of 
nonlinearity. This should be applied sparingly, as every time this rule is applied 
the gradient bounds deteriorate. It should definitely be applied at the 
initialization of the method. 
[0097] Even if the heuristic below is not applied, the B used should be upper 

triangular under some permutation of the variables. That is, for any constraint k 
and variable pair {ij} , either B HJ = 0 or B ¥ = 0 . This is to minimize the 

number of terms actually present in the quadratic expressions. 
[0098] To apply the rule in the initial step, we redistribute the symmetric 

elements in B according to the heuristic rule of putting all the quadratic "weight" 
on the variable with the smallest range (break ties lexicographically) — 
(Initialization at iteration 0.) 
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V{*,/,y}:if |v, - M .|<|v,.- M/ |,then j** ,+ V 

[0099] To apply the rule in subsequent step, we redistribute the symmetric 

elements in B only if there is a significant difference between the variables 1 
ranges — 

(At iteration m) 

V{k,iJ) : if |v. -u I < 10*|v, -iij and 5^ * 0, then T f = + ^ , 



Square-Free Bound Propagation 
[0 1 00] Assume we have a quadratic function inequality 

* /* to =Q+X 4«*i + Z - «** 

a set of point bounds w < x < v , and wish to refine the bounds for a variable x 3 . 
We will ignore any benefits of explicitly considering square terms x/ in the 
function f k , (hence "square-free propagation") although we will allow squared 
terms for variables other than Xj to appear. The Non-Square-Free Bound 
Propagation section explores the additional benefits of explicitly considering 

2 

Xj . 

[0101] In order to apply the Lemmas (see section below), we first rearrange the 

inequalities into the form 

lb k ~C k < A^Xj + B jjXj ) Xj +Yu A * x i + H B & x j x t 

and so if we define 



we also have 
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l K ~ C k ~ Z A * x > ~ Yu B m x ) x i ^r^ub k -C k - 24b*i " Z^*, x <- 

[0102] From the Bounding Lemmas section, we define the bounding functions 

for multiply and square — 

Vxy(x 0 ,x x ,y a ,y x ) = mmix^^xj^y^y,) 

ljxx{x a ,x x ) = min(x 0 2 , (x 0 x x ) + , x t 2 ) 
uu(x 0 , x, ) = max(x 0 2 , (x 0 x, ) + , x, 2 ) . 
We then define bounds fi 0 < j} < /?, and y 0 < y < y x for {fi,y} — 

Po = 4/ +Z(^ + *«,)♦«, + +^)_v f + (5,A"y + (5„)_v, 

A =4« v, +(*„)_«, 

r^^-Q-zaU-ZlAO.v, 

ivy 

- Z(%) + ^( M y' v y"«' v /)- Zfc)_Wfy(« y ,v y , Wl .,v,.) 
~ Z X . v, ) - Z ( B m )- » v « ) 

[0 1 03] Form the set {fi Q , # , p 0 , } . We now need to split up the set into the 

cases considered in Lemma B.4. 
[0104] Foreachset {/? 0 ,/?„y 0 ,^},ifit is the case that < 0 < fi x , then split 

the set into two sets, {fi 0 ,0, y 0 , ^ } and {0, # , y 0 , } . 
[0 1 05] For each set , # , ^ 0 , ^ } , if it is the case that y Q < 0 < ^ , then split 

the set into two sets, {/? 0 , ft , y 0 ,0} and {fi 0 , # ,0, ft } - 
[0 1 06] Up to four separate sets may result from this procedure. 
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[0 1 07] For each set B m - {ft 0 , # , x 0 , } > use Lemma B .4 to compute a range 

u f m < Xj < V m for each set. Restrict each range by the original bounds (allowing 
for the possibility that the range will be empty). 

u' m =max(w y ,0 
v' m =min(v J? vl) ' 

[0108] Given the resulting sets of ranges { {u f m , v' m } : m = 1 . M} , the variable will 

be in their union, x s e \J {u' m , V m } . 

m 

[0109] One can simplify them by sorting first by u f m then by v f m and then 

iteratively applying the following rule until no further simplifications are 
possible — 

1) If«;<«: +1 and < +1 <v;then 

{min(w^ M ; +1 Xmax(v^v; +1 )} = Kyju{u' m+l9 v' m+l } can replace both ranges. 

2) If v* m < u' m9 then delete the interval as being infeasible. 

[0110] Given a set of multiple simplified ranges, there are two possible ways to 

apply them to the original node — 

A) By conservatively bounding the multiple ranges, the node with bounds 
N = { i u j > v j } 9 {Uf j v J : i * J} can be refined to 

N = { {min m (u' m ),max m (v^ )} , {u t 9 v t }:i*J}. 

B) The node with bounds N = {{u J9 v J } 9 {u n v i }:i ' ± J) can be split into multiple 
refined nodes N m = {{u' m J m },{u n v t } :i*J). 

C) A third option would be to look at the benefit of splitting according to (B) 
over using the simple bounds of (A) by computing the ratio of the ranges 

!,(<-<) 



and only splitting the node if the benefit is sufficient, say p < .5 (the equivalent 
of one subdivision). 
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Non-Square-Free Bound Propagation 
[0111] The procedure here is similar to that in the Square-Free Bound 

Propagation section, except that the squared term x/ does appear in the 
function f k , so the defined bounds /? 0 < ft < p x and y 0 < y < y } will differ by a 
factor of the reciprocal of the coefficient of x/ , and Lemma B.5 will be invoked 
instead of Lemma B.4. 
[01 12] In order to apply the Lemmas, we first rearrange the inequalities into the 



form 



lb k -C k <A kJ x J + Ypu+BvK Xj+B ku x J 2 +Y d A H x i + <ub k -C k 

i*J i*J,j*J 



\i*J J 



and so if we define 

r 



4w 

V i*J J 



we also have 

T,B„,x J x t <B UJ r<ub k -C l - y £ t A u x l - Y B *XjXi 

[0113] We use the same bounding functions from the Square-Free Bound 

Propagation section. 

[01 14] We then define bounds B 0 < /? < /?, and y 0 < y < y x for {fi,y} based on 

the sign of B w — 

A = 4, + £(*„ +B w ) + u l +Y j (B kJi +B ldJ )_v i 
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- S ( B >* )+ > v . ) - 1] )- . v * ) 

*w i*j 

- X (** )+ » v , ) - X ( s k )- »«(». • v > ) 



5 U/ 



Yi =A- 



If5 uy <0 



A- A 



5 «/ 

A 



We now follow the same procedure as in the Square-Free Bound 
Propagation section, except that we apply Lemma B.5 to the generation of the 
ranges Xj e [JK, > v ^} . Note that unlike Lemma B.4, each case in Lemma B.5 
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may themselves generate multiple or empty ranges. 



Proofs for Nonlinear Function Problems 
[01 1 6] For the time being, let us consider only feasibility, although the same 

approach applies to optimization as well. 
[0117] Let f(x) : 9T -> be a function with Lipschitz-continuous 

derivatives. 

[0118] Assume we have a trial solution x . By the mean value theorem, for every 

point x there exists a point x (on the segment between x and x) such that 
(A.l) AW = A(x) + V/,(x)(x-x) 

[0119] Within a set of point bounds, 

(A.2) u < x < v 

there exist gradient bounds 

(A 3) Fk{UyV) = M * V/ * {X) :u ~ x ^ v ) 

G k (u 9 v) = sup{Vf k (x):u<x< v] 

so that 

(A4) u<x<v ^ F<Vf(x)<G 

(For notational convenience, we will suppress the functional dependence of F 
and Gonu and v.) 
[0 1 20] We will assume that the trial solution satisfies the bounds 

(A.5) u < x < v 

[0121] If we define 

u — u — x 
V = v-x 

(A.6) 

G=G~Vf(x) 
we will get the enveloping inequalities 
(A.7) u <x~x<v =^> 
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[0122] 



f{x) + Vf(x)(x -x) + F(x - x) + + G(x - x)_ 
<f(x)< 

fix) + V/(J)(x - x) + G{x - x) + + F(x - x)_ 
Now consider the following series of problems — 

f u°<x<v 0 1 



(A.8) P0 = jx 
(A.9) Pbd(u,v) = \x: 



lb<f{x)<ub\ 
u<x<v 



(A. 10) LP(x,u,v) = 



x,z,w: 



lb<f(x)<ubj 
u < x-x <v 

lb < f{x) + Vfix )ix -x) + Gz-Fw 
ub > fix) + Vfix)ix -x) + Fz-Gw 



x-x = z-w 

z + w< max(Jv|, \u\) 



z,w>0 

[0123] Define the infeasibility of a particular point to be 

(All) A k ix) = ma x(if k ix)-ub k ) + ,ilb k -f k ix)) + ) 
(As a result A(x) > 0, and A(x) = 0 only if the point x is feasible.) 
[0124] Theorem 1— 

Tl-1) For every xePO, there exist bounds {u,v) such that 

x<zPbdiu,v). 
Tl-2) Pbdiu,v)cLPix,u,v). 
Tl-3) For every {x,u, v} e LP(x,u, v) , 

A(x) < (G - F)(z + w) < (G - F) max(|4|v|) < (G - F)iv - u) . 
[0125] Corollary 1— 

Pbdiu,v)nLPix,u,v) = Pbdiu,v)nLPix',u,v), 
[0126] For the quadratic problem, 

(A. 12) fix) = C + Ax + Bxx 
and we can find 
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x =\(x + x) 
(A. 1 3) V/(x) = A=A + (B + B*)x 
W(x)-Vf(x) = B(x-x) 

so we can compute a set of gradient bounds 
G = (B) + v + (B)_u 

so that 

(A.15) M <x<v => F<Bx<G^> 
F = F-Bx < B(x-x) = Vf(x)-Vf(x) ^ G-Bx = G 

Bounding Lemmas 

[0127] Lemma B. 1— 

Let z = xy . If x 0 < x < x, and y 0 < y < y lt then 

mi «( x 0 y 0 >> 0 , ) < z < max(x Q y 0 , x, y 0 , x 0 y x , x l y x ) 
[0128] Lemma B.2— 

Let z = x 2 . If x 0 < x < x, , then 

min(x 0 2 ,(* 0 *i) + >*i 2 ) max(x 0 2 ,(x 0 x 1 ) + ,x I 2 ) 
N [0129] Or equivalents, Lemma B.3— 

Let z = x 2 . If x 0 < x < x, , then 

a) If x Q < x, < 0 or 0 < x 0 < x, then 
min(x 0 2 ,x! 2 ) < z < max(x 0 2 ,x, 2 ) 

b) If x 0 £ 0 < x t then 

0<z<max(x 0 2 ,x, 2 ) 
[0130] Lemma B.4— 

Let y=Bx.lf and y 0 < / ^ y } , then 

I)If 0<B 0 <B x and0<r o ^r, then 
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^-<x<^- 

A A 



(^<x<«>if J3 0 =0) 

Pi 



H) If fi 0 < P x <0 and 0 <y 0 < y„ then 



-^<x<-^- 
A Po 



(-™<x<^-ifp=0) 
Po 



HI) If 0<p o <p i and r 0 <^ <0 then 



-^<x<^- 
Po A 



(-co<*<^- if /? 0 =0) 



IV) If /? 0 < yS, <0 and ^ 0 <^ < 0, then 



(^-<x<oo if/?, =0) 



[0131] 



The following is an equivalent form, more convenient in some cases 
Lemma B.4' — 

Let y = p X . If p 0 < p < p i and y a <y<y„ then 
LIE) If 0<P o <P 1 then 



mm 



[Po'Pl) 



<x< max 



Yx Yx 



{Po'PxJ 



n,IV)If p 0 <p l <0, then 
< x < max 



mm 



r \ 
Yx Yx 



Po Px 



( Yo Yo" 



(±«coify? 0 =0) 

Po 



(i_«_ooif A=0) 

Pi 



A Ay 

[0132] Lemma B.5 — (Note that Cases I and II result in identical results.) 

Let x = ^x + x 2 .If p 0 <p<p l and y 0 ^y^y lt then 
I) If 0 < y? 0 < /?, and 0 < y 0 <y x then 



A \Px 



or more accurately, 



Po , Po 



+Yx 



either 
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4 



II) If A < (5 X < 0 and 0<^ 0 <^,, then 



_A_ A 2 ... A , /A 



4 " 2 V 4 



or more accurately, 



either 4-K +r , £I ,-A-K 



+r 0 



A 2 



A , A 



■+r, 



or + J— + r 0 <x<-^- + 

2 V 4 /0 2 , , 

111) If 0<j3 0 </S t and / 0 <x, <0 then 

Hl.a) If —j- + < 0 , then there are no possible solutions for x. 

B 2 

m.b)If ^- + y, >0,then 



2 V 4 
or more accurately, 



either -A- -j^+y, <x<-^--_ 



'El 
4 



+ 



IV) If A * A ^ 0 and y 0 < y x < 0 , then 
B 2 

HI.a) If + r , < 0 , then there are no possible solutions for x. 



0 2 

ELb)1f £ j- + y x >0,then 
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T~"VT 



or more accurately, 
either 



+r. 



or-« + 

2 ' 



+r 0 



Gradient Bound Subdivision 
[0133] It is also viable for the quadratic-only problem to subdivide the gradient 

bounds directly, e.g. by choosing a dimension {£',/'} and subdividing that 

F . + G . 

gradient bound, e.g. at// fe . = **' **' , or £T fe . = (Bx) h if subdivision at the 

linearization points is desired.. In that case, (1 .9) no longer applies, so it is 
necessary to include the gradient inequalities 

F <Bx<G directly in the problems (1 .23) - (1.28) as explicit inequalities. 
Theorem 1 and Corollary 1 will still apply. 
[0134] If one is subdividing the gradients, then there is additional bound 

propagation between the point bounds and the gradient bounds that can be 
defined — 

1) To refine the gradient bounds 

F' = max(F, (B)+u + (B)_v) 
G' = min(G, (B) + v + (B)_u) 

2) To refine the point upper bounds 



37 



AttyDktNo. H0002678/256.114US1 



V/:v;=v y 



Vk,i,j:< 



B ¥ >0^>v'j = min 



v' £*- 



B ¥ < 0 => Vj = min 



5, 



To refine the point lower bounds 

V/ : u } = ii . 



2^. > 0 = max 



j*j 



kji 



B^<Qz^ u'j = max 



Robust Linear Propagation 
[0135] The following modifies (2) in the Simple Propagation of Bounds to 

account for possible numerical problems. Let s > 0 be our numerical tolerance. 
For every linear constraint 

k:f k (x) = C k +J^A ki x i 

and for every variable x t which appears in f k and such that \A H \ > s , we define 

the following at each stage of the iteration where the bounds at that iteration are 
{u,v}. 

p»=i\-c k -Z<4M -S(4)-»; 



l + ZKl + max ( M ;|'K|)) 
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We can determine then determine new point bounds 



VkeK x , f k linear, V/ : 



( 



( 



A M > e => v' = max 



A M <-e^> v' = max 



Uj - l + ,min 



v., max 



Uj - [1 4- u } |, min v. , max 



J) 



VkeK x ,f k linear, Vi: 



■4b- > 5 : 



> ii ~ mm 



Vj 4- |l + v y |, max 



w ,mm 



V 



^ <-e^>u] -min 



v y + 1 + v^max 



K f? min 



Numerically Robust Square-Free Bound Propagation 
[0136] The following modifies the Square-Free Bound Propagation to account 

for possible numerical problems. Let e > 0 be our numerical tolerance. 
[0137] For convenience, define the following error bound on the product xy 

&xy(x 0 ,x l9 y 09 y l ) = e max(|x 0 1, \y 0 1, \x, |, \y, [) 

[0138] From the Bounding Lemmas, we define the robust bounding functions for 

multiply and square — 

II s xy(x 0 ,x 1 ,y Q9 y l ) = mm(x 0 y 0 9 x x y Q9 x 0 y x , x x y x ) - Axy(x 0 9 x X9 y 09 y x ) 
v € xy(x 0 ,x }9 y 09 y l ) = max(x o y 0 9 x x y 09 x 0 y x 9 x,y,) + Axy(x 0 9 x x 9 y 0 , y x ) 
fi £ xx{x 0 9 x x ) = min(x 0 2 , (x 0 x, )+,x x 2 )- Axy(x 0 ,x l9 x 09 x^ 
v e xx(x 0 9 x t ) = max(x 0 2 , (x 0 x x ) + ,x x 2 ) + Axy(x 0 9 x l9 x 09 x x ) 
We then define rigorous bounds /?* < J3 < fi* and y\<y< y[ for {fi>y}— 

A* = 4t/ + * + 2 < 5 « + ^ + £ > + < v * + 0 + 2 ( 5 « +B ul +e)_ («,. - £ ) + (Bjj + s) + (v, + s) - 
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Yl = lb k -C k - £ -Y l (A ki +e\( Vi + s)-Y J (A ki +s)_(u i -s) 



+s)yxx(u t ,v l )-Y J {B kii +e)_Li e xx(u n v l ) 



Yl = ub k -C k +5-2(4 "*)+(«, - £ )-Z(4 ~e)-(v,+e) 



We now follow the identical process as in (2.8) replacing the original set 

{A>>A>r 0 >r,}bythe set ipz,p; 9 rl,rl). 

Numerically Robust Non-Square-Free Bound Propagation 
[0139] The following modifies the Non-Square-Free Bound Propagation to 

account for possible numerical problems. Let s > 0 be our numerical tolerance. 

Define the robust bounding functions for multiply and square as in the 

Numerically Robust Square-Free Bound Propagation section. 
[0140] We define rigorous bounds J3* < J3 < fi[ and yl<y< y\ for 

{/?, y} based on the sign of B w — 

^^A kJ+£ + Y J {B Ui + + s) + (v, + e) + £ (2* w +B ku + e). (« f - e) 



i*J,j*J,i*j i*J,j*Jj*j 

~Y.i B m + *) + v ***("f> v /)-Zfoa +f)_^xx( M( .,v.) 
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ft =ub k -C k +e-% {A u - s\ (u, - s ) - £ (4. - s)_ (y t + s) 

[0141] We now follow the same process as in the Non-Square-Free Bound 

Propagation section replacing the original set {/?„ , /?, ,y 0 , y l } by the set 

Bound Propagation for Mixing Functions 
[0142] A common equation for refinery and chemical processes is a mixing 

equation of the form 

(G.l) x 0 (y x +y 2 +.. t + y n ) = Xiyi + Xi y 2 + .„ + x n y n 
where 

y x >0,y 2 >^..,y n >0. 
[0143] Here we are given materials indexed by {1,2, ...,«), and recipe sizes (or 

rates, or masses, or ratios, whatever mixing unit of interest) {y , , y 2 , . . . , y n ) for a 
mix of the materials. Then (G.l) represents a simple mixing model for a property 
x 0 of the resulting mix given values {x l9 x 29 ... 9 x m ) for the property of each of 
the materials. 

[0 144] In this case, the resulting property will be a convex linear combination of 

the individual properties. Hence if we are given bounds 
(G.2) a, <x x <b v a 2 <x 2 <b 2 ,„.,a n <x n <b n 
we can derive the bounds 

(G3) X ° ^ min ( a i>*2>».>tf«) 
x^max^,^,,..,^)' 

[0145] It is to be understood that the above description it is intended to be 

illustrative, and not restrictive. Many other embodiments are possible and some 
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will be apparent to those skilled in the art, upon reviewing the above description. 
For example the local linear bounding and the local linearization can be 
performed in the opposite order, and more. Therefore, the spirit and scope of the 
appended claims should not be limited to the above description. The scope of 
the invention should be determined with reference to the appended claims, along 
with the full scope of equivalents to which such claims are entitled. 
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